# dvir, 2011-08-30
# this class is for cases when we need to create many dat files, by a given dat file.
# the dat file will be read only once to the memory !

# this class isn't perfect for all dat's, but it can be modified.
# for now, I'm sure that it's works with : sex-linked:no, no liability classes, no complex loci (=not a disease of two loci)
package DatWrapper;

# details about the structure of a dat file can be found at : 
# http://bioinfo.cs.technion.ac.il/superlink/Locus_File.html

use strict;
use General;
use Distances;

use Class::Struct;

my $ONE_LOCI_DISEASE_CODE = 0;
my $TWO_LOCI_DISEASE_CODE = 2;

my $REG_MARKER_CODE = 3;
my $DISEASE_MARKER_CODE = 1;
my $COMPLEX_MARKER_CODE = 4; # works with $TWO_LOCI_DISEASE_CODE

struct MarkerData =>
[
	name => '$',
	geneFrequencies => '@', # the gene frequencies of the marker
	distIn_Theta_FromPrevMarker => '$',
	distIn_cM_FromPrevMarker => '$',
	distIn_cM_FromFirstMarker => '$',	
];

struct PenetranceClassData =>
[
	posibilityForAffectionWhen_H_H => '$',
	posibilityForAffectionWhen_H_D => '$',
	posibilityForAffectionWhen_D_D => '$',
];

sub new( $ ) {
	my $class = $_[0]; # the name of this class
	# the dat file seperator can be tab or space or both.
	my $datFileName = $_[1]; # may include a path
	
	my %self = (); # empty hash for the class
		
	open DAT_FILE, $datFileName;
	# line1
	# 36 0 0 8 0 << num of markers, ignored, sexLinked, superlink program, one/two locus disease
	# num of markers -> 1 more than the total num of markers, cause it includes the disease marker
	# one/two locus disease : 
	#	2 -> 2 loci disease
	#	0 -> 1 loci disease
	my $line1 = <DAT_FILE>;
	
	my @line1Arr = split /\s+/, $line1;
	my $numOfMarkers = $line1Arr[0]; # including the diesase marker
	$self{IS_SEX_LINKED} = $line1Arr[2];
	$self{PROGRM_CODE} = $line1Arr[3];
	$self{NUM_OF_DISEASE_LOCI} = $line1Arr[4];
	
	# making sure that this is only one marker,
	# sihce for now the system supports only one disease locus
	if ($self{NUM_OF_DISEASE_LOCI} != $ONE_LOCI_DISEASE_CODE) {
		#for omer request, we now support the possibility of two marker diesase
		#&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "We doesn't support dat file of more than 1 loci disease");
	}
	
	# line2
	<DAT_FILE>; # 0 0 0 0 << ignored
	# line3
	<DAT_FILE>; # the line is a list of markers numbers
	
	# reading the markers allels data
	my $currMarkerNumWithDieseaeMarkers = 1;
	my $currMarkerNumWith_Out_DieseaeMarkers = 1;
	while (my $currLine = <DAT_FILE>) {
		my @currLineArr = split /\s+/, $currLine;
		my $markerType = $currLineArr[0];		
		
		if ($markerType == $REG_MARKER_CODE) {
			
			# checking that the marker has a name
			if ($currLine !~ /#(.+)#/) {
				&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "marker no' $currMarkerNumWith_Out_DieseaeMarkers, doesn't has a name.\n$currLine");
			}
			my $markerName = $1;
		
			my $numOfAllels = $currLineArr[1];
			my $geneFrequenciesLine = <DAT_FILE>;			
			my @geneFrequenciesArr = split /\s+/, $geneFrequenciesLine;
			# making sure that we have enough frequencies
			for (my $currFreq = 0; $currFreq < $numOfAllels; ++$currFreq) {
				# checking that we have a frequency as a number
				if ($geneFrequenciesArr[$currFreq] !~ /\d+|\d+\.\d*/) {
					&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "Wrong num of allel frequncies for marker no' $currMarkerNumWith_Out_DieseaeMarkers.\n$currLine");
				}
			}
			
			# checking that we don't have too many frequncies
			if (defined $geneFrequenciesArr[$numOfAllels] and $geneFrequenciesArr[$numOfAllels] =~ /\d+|\d+\.\d*/) {
				&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "Wrong num of allel frequncies for marker no' $currMarkerNumWith_Out_DieseaeMarkers (Too much frequncies).\n$currLine");
			}
				
			# updating the markers data structure
			my $markerData = MarkerData->new();
			$markerData->name( $markerName );
			@{$markerData->geneFrequencies} = @geneFrequenciesArr[0..($numOfAllels-1)];
			
			# each entry in the array will be MarkerData, the array is from cell no' 1 and not from cell no' 0.
			$self{MARKERS}[$currMarkerNumWith_Out_DieseaeMarkers] = $markerData;
			++$currMarkerNumWith_Out_DieseaeMarkers;
		}
		elsif ($markerType == $DISEASE_MARKER_CODE) {
			# reading the next line
			$currLine = <DAT_FILE>;
			@currLineArr = split /\s+/, $currLine;
			$self{NORMAL_GENE_FREQUENCY} = $currLineArr[0];
			$self{MUTANT_GENE_FREQUENCY} = $currLineArr[1];
			
			$currLine = <DAT_FILE>;
			@currLineArr = split /\s+/, $currLine;
			$self{NUM_OF_LIABILITY_CLASSES} = $currLineArr[0];
			
			# in case of sexLinked chromosome, then the num of peneterance classes is doubled,
			# since that each sex has its one penetrance class
			my $numOfLiabilityClasses = $self{NUM_OF_LIABILITY_CLASSES};
			if ($self{IS_SEX_LINKED}) {
				$numOfLiabilityClasses = 2 * $self{NUM_OF_LIABILITY_CLASSES};
			}
			
			for (my $currPenetranceClass=0; $currPenetranceClass<$numOfLiabilityClasses;++$currPenetranceClass) {
				$currLine = <DAT_FILE>;
				@currLineArr = split /\s+/, $currLine;
				my $currPenetranceClassData = PenetranceClassData->new();
				$currPenetranceClassData->posibilityForAffectionWhen_H_H ($currLineArr[0]);
				$currPenetranceClassData->posibilityForAffectionWhen_H_D ($currLineArr[1]);
				$currPenetranceClassData->posibilityForAffectionWhen_D_D ($currLineArr[2]);
	
				$self{PENETRANCE_CLASSES}[$currPenetranceClass] = $currPenetranceClassData;
			}
		}
		elsif ($markerType == $DISEASE_MARKER_CODE) {	 
			&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "not supported right now");
		}
	
		# exiting after reading all the markers data
		if ($currMarkerNumWithDieseaeMarkers == $numOfMarkers) {
			last;
		}
		++$currMarkerNumWithDieseaeMarkers;
	}
	
	# reading the recombination information
	<DAT_FILE>; # this line is the: sex difference,  interference 
	my $recValuesLine = <DAT_FILE>;	
	
	#Omer 15/9/2011 - reading the program-specific options
	my $programOptions = '';
	while (my $currLine = <DAT_FILE>) {$programOptions .= $currLine;}
	$self{PROGRAM_OPTIONS} = $programOptions;
	close DAT_FILE;
	
	#Compute marker positions
	my @distancesInThetaForEachMarkerFromItsPreviousMarkerArr = ();
	my $numOfMarkersWithoutDiesaseMarkers = $numOfMarkers - 1;
	&Distances::calcDistInTehtaForEachMarkerFromItsPreviousMarker(\$recValuesLine, $numOfMarkersWithoutDiesaseMarkers, \@distancesInThetaForEachMarkerFromItsPreviousMarkerArr);
	
	#Save marker positions inside the object
	my $currDistIn_cM_FromFirstMarker = 0;
	my $currMarkerNum = 1;
	for my $currDist (@distancesInThetaForEachMarkerFromItsPreviousMarkerArr) {
		my $currDistIn_cM = &Distances::recToCm($currDist);
		$currDistIn_cM_FromFirstMarker += $currDistIn_cM;
		$self{MARKERS}[$currMarkerNum]->distIn_Theta_FromPrevMarker($currDist);
		$self{MARKERS}[$currMarkerNum]->distIn_cM_FromPrevMarker($currDistIn_cM);
		$self{MARKERS}[$currMarkerNum]->distIn_cM_FromFirstMarker($currDistIn_cM_FromFirstMarker);
		++$currMarkerNum;
	}
		
	return bless \%self, $class;
}

# since some executeable (for example haplo lars) requires a specific program code
# than we are changing the program code.
# Omer 15/9/2011 - the second parameter is the (optional) program-specific options
sub setProgramCode( $ $ ) {
	my $class = $_[0]; # the name of this class
	my $givenProgramCode = $_[1];	
	$class->{PROGRM_CODE} = $givenProgramCode;	
	$class->{PROGRAM_OPTIONS} = (defined($_[2])) ? $_[2] : '';
}

# writes to the given file handle, a dat file only of the given list of markers.
sub writeSubDatFileByListOfMarker( \* \@ ) {
	my $class = $_[0]; # the name of this class
	# reference for file handle to which the dat file will be written (concat at the end of the file)
	my $destFileHandleRef = $_[1];
	# assuming that the list of markers is a sorted lisr
	my $listOfMarkersArrRef = $_[2];

	my $numOfMarkers = scalar(@$listOfMarkersArrRef);
	my $numOfMarkersPlus1 = $numOfMarkers + 1;
	
	my $newLine1 = "$numOfMarkersPlus1 0 $class->{IS_SEX_LINKED} $class->{PROGRM_CODE} $class->{NUM_OF_DISEASE_LOCI} << num of markers, ignored, sexLinked, superlink program, one/two locus disease\n";
	my $markersNumbers = join ' ', 1..$numOfMarkersPlus1;
	print $destFileHandleRef 
		$newLine1. # line1
		"0 0 0 0 << ignored\n". # line2
		"$markersNumbers << markers order\n". # line3
		"$DISEASE_MARKER_CODE 2 << disease marker\n". # line4		
		"$class->{NORMAL_GENE_FREQUENCY} $class->{MUTANT_GENE_FREQUENCY} << disease frequency (normal gene freq, mutant gene freq) \n". # line5
		"$class->{NUM_OF_LIABILITY_CLASSES} << number of liability classes\n";


	# printing the penetrance class data
	for my $penetranceClass (@{$class->{PENETRANCE_CLASSES}}) {
		print $destFileHandleRef 
			$penetranceClass->posibilityForAffectionWhen_H_H().' '.
			$penetranceClass->posibilityForAffectionWhen_H_D().' '.
			$penetranceClass->posibilityForAffectionWhen_D_D()." << penetrance classes\n";
	}
	
	# printing the markers data & checking that the markers are sorted and all the markers exists
	my $prevMarkerNum = 0; # for start
	# we reduce 1 from the num of markers, because cell 0 in the markers array is empty.
	my $numOfMarkersInOrigDat = scalar (@{$class->{MARKERS}}) - 1;
	for my $currMarkerNum (@$listOfMarkersArrRef) {
		# checks
		if ($currMarkerNum > $numOfMarkersInOrigDat || $currMarkerNum <= 0) {
			&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "wrong marker num: $currMarkerNum, numOfMarkersInOrigDat: $numOfMarkersInOrigDat");
		}
		elsif ($currMarkerNum <= $prevMarkerNum) {
			&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "The markers list isn't sorted!\ncurrMarkerNum: $currMarkerNum, prevMarkerNum: $prevMarkerNum");
		}
		$prevMarkerNum = $currMarkerNum;
		
		#prints
		my $name = $class->{MARKERS}[$currMarkerNum]->name();
		my @geneFrequencies = @{$class->{MARKERS}[$currMarkerNum]->geneFrequencies()};
		my $numOfGeneFrequncies = scalar(@geneFrequencies);
		
		print $destFileHandleRef
		 "$REG_MARKER_CODE $numOfGeneFrequncies #$name# << allele numbers, no. of alleles\n".
		 join(' ', @geneFrequencies)." << marker frequency\n";
	}
	
	# printing the recombination values data
	print $destFileHandleRef 
		"0 0 << sex difference, interference\n".
		"0"; # the first recombination value
		
	# passing over all the cells, except for the first cell (cell 0) (cause its distance is always : 0)
	# the last cell index is ($numOfMarkers-1)
	for (my $currMarkerIndexInList=1;$currMarkerIndexInList< $numOfMarkers; ++$currMarkerIndexInList) {
		my $dist = $class->calcDistIn_Theta_BetweenTwoMarkers(
			$$listOfMarkersArrRef[$currMarkerIndexInList-1],
			$$listOfMarkersArrRef[$currMarkerIndexInList]);
		
		print $destFileHandleRef " $dist";
	}		
	print $destFileHandleRef " << marker distances\n";
	
	# printing the program-specific options
	if ($class->{PROGRM_CODE} == 8) {
		print $destFileHandleRef "1 << number of disease loci\n";
		print $destFileHandleRef
			"-n 2 1 $numOfMarkers << REC VARIED, INCREMENT, FINISHING VALUE\n";
	}
	else {print $destFileHandleRef $class->{PROGRAM_OPTIONS};}
}

# computes the distance in Theta between 2 markers
sub calcDistIn_Theta_BetweenTwoMarkers() {
	my $class = $_[0]; # the name of this class
	my $firstMarkerIndex = $_[1];
	my $secondMarkerIndex = $_[2];
	
	my $numOfMarkers = $class->getNumOfMarkers();
	if ($firstMarkerIndex <= 0 || $secondMarkerIndex <= 0 || $firstMarkerIndex > $numOfMarkers || $secondMarkerIndex > $numOfMarkers) {
		&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "calcDistInTehtaBetweenMarkers: bad indexes: firstMarkerIndex: $firstMarkerIndex, secondMarkerIndex: $secondMarkerIndex, numOfMarkers: $numOfMarkers.");
	}
	
	my $distIn_cM_BetweenTwoMarkers = $class->calcDistIn_cM_BetweenTwoMarkers($firstMarkerIndex, $secondMarkerIndex);

	my $distIn_Theta_BetweenTwoMarkers = &Distances::cmToRec($distIn_cM_BetweenTwoMarkers);
	return sprintf( "%.16f", $distIn_Theta_BetweenTwoMarkers);
}

# computes the distance in cM (=centi morgan) between 2 markers
sub calcDistIn_cM_BetweenTwoMarkers() {
	my $class = $_[0]; # the name of this class
	my $firstMarkerIndex = $_[1];
	my $secondMarkerIndex = $_[2];
	
	my $numOfMarkers = $class->getNumOfMarkers();
	if ($firstMarkerIndex <= 0 || $secondMarkerIndex <= 0 || $firstMarkerIndex > $numOfMarkers || $secondMarkerIndex > $numOfMarkers) {
		&General::reportOnErrorReleaseSemaphoresAndDie("DatWrapper.pm", "", "calcDistIn_cM_BetweenTwoMarkers: bad indexes: firstMarkerIndex: $firstMarkerIndex, secondMarkerIndex: $secondMarkerIndex, numOfMarkers: $numOfMarkers.");
	}
	
	my $distIn_cM_BetweenTwoMarkers = 
		$class->{MARKERS}[$secondMarkerIndex]->distIn_cM_FromFirstMarker() - 
		$class->{MARKERS}[$firstMarkerIndex]->distIn_cM_FromFirstMarker();

	return $distIn_cM_BetweenTwoMarkers;
}

# returning how many markers are in the given dat file
sub getNumOfMarkers() {
	my $class = $_[0]; # the name of this class
	# because cell 0 in the markers array is empty, than we are reducing 1.
	my $numOfMarkers = scalar (@{$class->{MARKERS}}) - 1;
	
	return $numOfMarkers;
}

# returning the marker name
sub getMarkerName() {
	my $class = $_[0]; # the name of this class
	my $markerNum = $_[1];
	
	if (defined $class->{MARKERS}[$markerNum]){
	return $class->{MARKERS}[$markerNum]->name();
}
	print "$markerNum is nor defined\n";
	return "$markerNum is nor defined";
}

sub getDistIn_cM_FromPrevMarker() {
	my $class = $_[0]; # the name of this class
	my $markerNum = $_[1];
	
	return $class->{MARKERS}[$markerNum]->distIn_cM_FromPrevMarker();
}

sub getDistIn_cM_FromFirstMarker() {
	my $class = $_[0]; # the name of this class
	my $markerNum = $_[1];
	
	if (defined $class->{MARKERS}[$markerNum]){
	return $class->{MARKERS}[$markerNum]->distIn_cM_FromFirstMarker();
}
	print "$markerNum is nor defined\n";
	return "$markerNum is nor defined\n";
}

1;